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Abstract: The XY-model shows in two dimensions in the strong coupling regime 
a universal distribution, named BHP, which in turn also describes other models 
of criticality and self-organized criticality and even describes natural data as river 
level and flow. We start by analysing the two dimensional XY-model and calculate 
the BHP probability density function. The results obtained for several dissimilar 
phcnomenons which includes the deseasonaliscd Danube height data raised the uni- 
versality hypothesis for rivers. This hypothesis is tested for the Iberian river Douro. 
Deviations from the BHP are found especially for medium and small runoffs. For 
regimes closer to the natural flow the fluctuations tend to follow the universal curve 
again. 

PACS numbers: 05.40.a, 05.65+b, 64.70.Nd 



1 Introduction 

In their pioneer paper the authors Bak, Tang and Wiesenfeld, [10], proposed 
the hypothesis that under very general conditions nonequilibrium systems 
consisting of many interacting constituents may exhibit universal behaviour. 
This behaviour among other features is characterised by the formation of 
correlations and should occur naturally without any parameter tuning from 
the outside. For these systems the dynamical response is complex but the 
statistical properties are governed by simple power laws. Since only dimen- 
sionless numbers can be raised to arbitrary powers in a meaningful way, 
this self-organisation of the system implies independence of any particular 
scale of length or time. Hence the behaviour of the systems on these self- 
organized states is close to those of equilibrium systems at the critical point. 
This phenomenon is known as self-organized criticality (SOC). 
In equilibrium thermodynamics criticality is linked with continuous phase 
transitions, [19]. The phase transition is continuous when the rate of change 
of some order parameter, and not the order parameter itself, is discontinu- 
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ous at the critical point. In these systems when temperature is equal to the 
transition temperature the spin correlation function instead of an exponen- 
tial decay reveals a power law behaviour and the corresponding exponent 
is called the critical exponent, in this case for the correlation function. As 
a consequence any local perturbation can be propagated through out the 
entire system and hence any member of the system affects all the others. 
At the critical point the contribution of the interaction between widely sep- 
arated points for the large scale fluctuations of the order parameter are not 
exponentially rare. In this frame the central limit theorem is not applicable 
and Gaussianity is not present at all. 

The existence of power laws statistics and critical exponents raises the ques- 
tion of classification of these systems. The phenomenon, whereby dissimilar 
systems exhibit the same critical exponents, is called universality. Two 
systems are assigned to the same universality class if they share the same 
dimensionality d of the underlying lattice (number of spatial dimensions) 
and the same dimensionality of the order parameter, D. All the systems in 
the same universality class have the same critical exponents. 
In this work we follow Bramwell, Holdsworth, Pinton, [H] in calculating the 
probability density function (PDF) for the fluctuations of the magnetic or- 
der parameter of a two-dimensional model (2dXY) for spins in the strong 
coupling (low temperature) regime using the spin wave approximation where 
it shows a universal distribution, i.e., independent of system size and crit- 
ical exponent v. This universal PDF describes other models of criticality 
and SOC [7] hence we show the data collapse for the deseasonalized Danube 
river height data at Nagymaros (80 years) following the steps of [3j and com- 
pare it with the daily mean runoff data of river Douro at Regua (26 years). 
According to those authors no important differences occur when height or 
runoff measures are used. 

This article is organised as follows the second section describes the two- 
dimensional Ising model (2dlsing) and its frame set as the preliminary for 
the 2dXY calculations presenting, order parameter and PDF's at different 
regimes of temperature; subcritical critical and supercritical. In the third 
section we present the analytics for the spin-wave (SW) approximation of 
the 2dXY model and the magnetic order parameter PDF for the critical (low 
temperature) regime named BHP, after Bramwell, Holtsworth and Pinton 
[7]. In the fourth section we show the data collapse of the deseasonalized 
Danube river height data and we compare it with the Douro river in two 
cases, the entire year and the winter regime. Apparently the deviations of 
the Douro river from the BHP are due to some river flow regulations in 
order to have fluctuations closer to the Gaussian PDF. For larger values of 
the streamflow the regulation is no longer possible due to storage limitations 
and the river dams are set open. For larger values the river fluctuations are 
much closer to the BHP form. The main differences between the Danube 
and Douro are related to the amount of water on the river basin and its 
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distribution in time. The Douro is a southern European river where the 
summer is usually very dry contrasting with the Danube basin where there 
is lots of water the year around and regulation seems not to have any major 
effect on the natural flow. 



2 The Ising model 

Ernst Ising suggested a very simple thermodynamic model to understand 
spontaneous magnetization (Ising, 1925, [K]). Statistical mechanics as used 
by Ising considers equilibrium distributions 

p*(a 1 ,...,a N ) = ±---e H ^°°>->^ (1) 

and quantities derived from them, like the magnetization per spin (M) := 
(lEiIi (I «) an d magnetic susceptibility x := (M 2 ) — (M) 2 . Here the N 
spins {(Tj}^, Oi 6 { — 1,+1}, are distributed in p*{<7\, erjv) according 
to Boltzmann weights e w (originally introduced to represent well known 
distributions in ideal gas theory, the Maxwell distribution of velocities) with 
a function 

N N N N 

H(ai, a N ) = ^2 E V ij a i a j + E hm + E Ci ( 2 ) 

i=l j = l i=l i=l 

and a normalization, the partition function Z. In order to obtain reasonably 
sized numbers we use C = J2iLi c * = ~~ Nha.2 for the 2^ possible configura- 
tions of spins up and down, hence 

Z= ^] ... e^= 1 ^i= lVijai(7j+ ^= lhiai ■ . (3) 

IT1=±1 (Tjv=il 

It is Vij := V ■ Jij with coupling strength V and adjacency matrix J^- G 
{0, 1} and h := hi the external magnetic field. Then the distribution of the 
magnetization per spin is given by 

p(M)= £ ... « M-i^ ffi ]/K... )ffff ) (4) 

CT1=±1 CTJV=±1 V 1=1 / 

from simply applying Bayes' rulep(M) = £^=±1 ••• E CTAr =±i^( M l cr i' • •> CJ A r ) 
p(cri, ...,<tjv) for joint and conditional probabilities. 

The mean magnetization, (M) = Y^L Q M-p(M), with M = -1+i-AM, 
AM = 2/N is given by 

(m}= - E (E M -M M -^E^))^'-'^) • ( 5 ) 

<ti=±1 <tjv=±1 \i=0 V i=l / / 
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Using the <5-function the mean magnetization is then 

<">- E •• E fef:*)^ • < 6 > 

<T1=±1 CTjv=±l \ 1=1 / 

The mean magnetization (M) is an easily accessible quantity in physical sys- 
tems, whereas the distribution of the magnetization p(M) is easily accessible 
in small systems. Also the mean (M) is dominated by the maximum M max 
of the distribution p{M) for large systems. The Ising model was proposed 
by Lenz|15j and solved by Ising for one dimension. The two dimensional 
case shows already a phase transition. In 1944 Onsager|17|. solved the Ising 
system for d = 2 in the absence of an external magnetic field. Here we 
present a numerical simulation for the magnetization distribution for the 2d 
Ising model for a N = 4 x 4 lattice with different coupling strength, V. The 
external magnetic field is absent, h = 0. 




Figure 1: Coupling 
strenght value is V — 0.13. 
The system is subcritical, 
only one maximum at M = 
appears. 
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Figure 2: For the value 
of V = 0.14 « V c we 
find one maximum at M = 
0, but being broadened up 
just before splitting into 
the two symmetric maxima 
M max / 0. 
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Figure 3: Coupling 
strenght value is V = 
0.145 two maxima Mmax 7^ 
are found symmetrically 
around M = 0, where 
there is a local minimum in 
p(M). 



For low temperature (strong coupling), fig. [T] the spins are strongly cor- 
related and tend to align with their neighbours so the distibution mode is 
for zero magnetization. Near critical temperature, fig. [2] small variations 
of spins can be propagated through the lattice and will affect spins at large 
distances. That is the point of phase transition. Fluctuations can easily 
happen for values between M ~ —0.5 and M ~ +0.5, for which all magne- 
tizations in this range are about equally likely. For temperatures just above 
criticality, fig. O two maxima start appearing symmetrically around the 
origin. 
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2.1 The XY- model in two dimensions 

In the XY model spins are able to rotate freely in the XY-plane, i.e. the 
spin variables are with 

where the angle 9i changes between —it and +n. 
The stationary distribution is again given by 

P*(S 1 ,...,S N ) = ±-e^-^) ( 8 ) 

and 

N N N 

ns,, ...,s N ) = E E + E^ + c ( 9 ) 

i=i j=i i=i 

where later the constant C will be fixed to obtain reasonably small values 
for the partition function Z (irrelevant for all physically relevant quantities, 
e.g. the stationary distribution p*). 

Since the spins are unit vectors rotating we can replace them just by 
angles between two interacting spins 

SiSj = \Si\ ■ \Sj\ ■ cos(9 l - dj) = cos{9i - 9j) (10) 

or the angle between spin and outer magnetic field (see |19} I22j) 

HS i = \H\-\S i \-cos{6i) = h-cos(e i ) . (11) 

Hence 

p *(e x ,...,e N ) = ^-e H ^>~M ( i2) 

and 

JV JV N 

H(9 1 ,..£ N ) = vJ2zZ J y cos ^ i ~ 9 ^ + zZ h ' ms ^ + c ( 13 ) 

i=l j=l i=l 

with partition function as normalization constant 

Z= dBt... r d9 N e n ^-M . ( 14 ) 

J —71 J —TV 

The absolute value of the magnetization as order parameter for the phase 
transition [18] is given by 



m 



1 N 1 

_L ii v 9 n= — 



N 



E cos (^ 



+ 



vi=l 



' N 

E 

\i=i 



sin[ 



--: M(9 l: . 



> On) 
(15) 
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and the distribution of this quantity 

p(m)= d0i... / dffjv ^(m-Af^i,...,^)) -e w ^--^ . (16) 



This defines completely the XY-model in any dimension (where the dimen- 
sion is fixed by specifying the adjacency matrix J in On)). This will 
be referred to the non-quadratic model, since later for analytical treatabil- 
ity of the strong coupling regime the cos-interaction will be replaced by the 
quadratic Taylor expansion (valid for small angles between the spins only or 
spins and external magnetic field). The BHP-distribution will turn up by 
evaluating the Eq. (|16|) in the case of quadratic approximation of the above 
model. Using the second order Taylor series expansion for the cosine, 

cos(x) = 1 - -^x 2 + 0{x 4 ) (17) 

hence we can replace in Tt(0i, #/v) in the strong coupling case (V » 1) 

cosiOi-O^^l-^Oi-Oj) 2 (18) 

i.e. approximating all relevant quantities in quadratic form in the exponen- 
tial function and then use Gaussian integration. This approximation is valid 
for strong coupling V >> 1, where in the physical model spin waves are the 
only relevant dynamical phenomena. 
The partition function as normalization constant 

Z = I dOi... I dO N e H(di,-fi N ) ( 19 ) 

J —TT J —TV 

becomes in quadratic approximation 

/OO /'OO f"K 

d0 1 ... / dO N ^... / dO N e J *°*- v e - v - e - trM - (20) 
-OO J —OO J —TT 

It is A := Q ■ 1 — J with J the adjacency matrix, Q := Y^f=i Jij ^ ne number 

of neighbours, and Jtot '■= Sj=i Ylf=i Jij * ne total number of connections in 
the adjacency matrix. Changing variables y = TO with the transformation 
matrix 

T := (u nk ) = (-Le- 2 ^A (21) 



diagonalizes matrix A such that AT = A with A the diagonal matrix of 
eigenvalues of A and u k the eigenvectors of A due to the structure of the 
adjacency matrix J. T^ is the Hermitean matrix of the complex matrix T, 
i.e. the transposed and complex conjugate. In 2 dimensions 
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with k x etc. having values 1,2, ..,L, and with N = L x L, the variable 

k:=k x + {k y — 1) • L (23) 
has values 1,2,..., N. Likewise for n. Also in 2d the eigenvalues are 



A fc = 4 - 2 • cos ( 2-n—kx j - 2 • cos \2ir-ky ) . (24) 



In the general case the partition function is then given in changed coordi- 
nates, the eigencoordinates of A, 

/oo poo rc 

d Vl ... dy N -!... dy N e- v -y trA y-\det(T)\ (25) 
-oo J —oo J —c 

with \det(T)\ = 1 and transformed integration boundaries for the zero eigen- 
value coordinates c. It turns out that the N th coordinate is the zero eigen- 
value coordinate, in which case c = (wjvAt) -1 ' 71 " = • 7T - 
Performing the Gaussian integratioro of Eq. gives 

N-i 



z = eJtot ' v '\\ \hrrv ' ^ 



k=l 1 K 



where the last factor comes from the eigenvalue zero. Using the Fourier 
representation of the delta function in [161 we obtain 



r*oo pn poo -i -i 

P (m)= I ■■■ I / de x ---de N ± e i(m-M(e)x dx ^ e H(e 1 ,-,e N ) _ 



2-7T 

— OO J — TT J —OO Aj " 



(27) 

For simplicity purposes it's more convenient to write the magnetization in 
terms of the spin variables instead of vector spins, 



M = iE^- ) 2 ( 28 ) 

r=l 

this is just the average, (6 r ) for unconfined spins. For the case of periodic 
spins, following [9], the instantaneous magnetization direction, 8, is defined 
by, 

= arc tan ( E^l^L ) ( 29 ) 
\Er=l cos r J 



Let B be a complex symmetric matrix with a non-negative real part and non zero 
eigenvalues, then, f °° ■ ■ ■ f°° d6 e~- B - - 

° ' J— oo J— oo — 
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which is the same as the average for large N. Then we may write M(9) = 
1 — 2]v#* r #. The magnetization PDF can be written using the quadratic 
approximation , 



p(m ) = _±1 I"" g f°° . . . f°° f dQ^-i+^eyve-^ (30) 



r A9 



where d9 = Ylk=i d@k- Simplifying for Gaussian integration, 



Z J — OO J - X ./ - x ' - T 



Using the linear transformation (|2ip and writing i? = V.A — 4^ 1 there is a 
diagonal matrix A such that B = TKT^. Using the multivariate Gaussian 
integral, we obtain 



/oo poo / N 

The matrix B is equal to q(A) where q() is a polynomial of A. Then we may 
write, 

N N 



det(B) = Hl k = ~[[q(\ k ) (33) 
k=l k=l 

eigenvalues of hie 
values of B. The magnetization PDF is, 



where {\k}k=i are the eigenvalues of matrix A and {(-k\k=i are * ne e ig en_ 



(34) 



where c = vNtt. Simplifying we have, 

(N—l \ — 

n • <-) 

The mean and standard deviation magnetization can be calculated using 
the distribution of the magnetization absolute value, [3H Under the same 
assumption for the A matrix we have 



(to) — 1 ^ Y,k=l 2VX k 
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At this point one has to do the appropriate normalisation in order to get the 
universal critical PDF. The distribution function for the magnetization has 
been discussed by several authors (see references in [9]). Binder, 1992|20j. 
argued that from the Ising model with distribution f(m) there is scaling of 
the form 

p(m) ~ L^^pimL 13 ^ , L/C) (37) 

where £ is the correlation length. We know that very close to the critical 
point the correlation length is going to be much bigger than the system size. 
In this case L is the only important scale so / should depend only on the 
variable mL^l v , 

p(m) ~ #p(m#) . (38) 

This form is not surprising since that for this system the mean (m) scales 
with L~P/ U and being critical the standard deviation, a is also scales in 
the same way as the mean value. These facts explain the universal form of 
p(mLPI v ). Hence we can rewrite equation (|38p as 

p(m) -—/(—) • (39) 

The standard deviation is the correct normalization for the order parameter. 
Defining [i as 

fj, = (40) 

it is straightforward to see that the PDF p(fi) = a m p{m) and this PDF is 
normalized not depending on system size 

a m .p(m)~L Q .f(^-^.L°) . (41) 
The form of the normalized and consequently universal BHP is, 



00 dx 

-oo 27r 



1 ^ 1 W^EtY^-Efji 1 [^^-iarctan( 

2AT2 2^ A 2 e 
k=l k 



.e 



lnf 



1+ 



k 



(42) 



In figured] we show a BHP PDF simulation for L = 10, iV = L 2 , using eq. 

(021). 

At this point we should mention that the magnetic order parameter density 
of the 2dXY model approaches this universal PDF when T — > 0. Surpris- 
ingly the quadratic approximation remains close to the true PDF for almost 
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Figure 4: Logarithm of the universal BHP PDF for the 2dXY model with L = 10. 



all the critical regime. It has been shown recently |14[ [TT] that the weak 
temperature dependence gets noticeable near the KT phase transition. 

3 BHP and the Danube and Douro data 

Janosi and Gallas [2 J analysed statistics of the daily water level fluctuations 
of the Danube collect over the period 1901-97 at Nagymaros, Hungary. The 
authors found in the Danube data similar characteristics to those of company 
growth [15]. This suggested that a universal description of the statistics 
should exist for both systems. Those authors noticed a data collapse for the 
conditional PDF of the one day logarithmic rate of changes for the mean 
adjusted river height. 

Following [2J, let h(t) be the 87x365 data points^ the time series of height 
measurements and define a season mean and, following [6], a season standard 
deviation, h(t) and ah, respectively, 




(43) 



3=0 



86 
3=0 



-i 365 



where h(t,j) = {h(t + j * 365)} 



J t=i 



2 th 
all the 29 of February were discarded 
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^M = ^kpZ7^f . (44 ) 

The seasonality adjustment consists in normalizing each value of the original 
time series by subtracting the respective component of the season mean 
and dividing by the respective component of the season standard deviation, 
(^r) ' wa t er plays for the river the same role as the coupling strength 
in the 2dXY model. In figures 06] we show the Histogram and Log histogram 
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Figure 5: Histogram of the seasonally 
adjusted daily water height Danube river 
data. 
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Figure 6: Log histogram of the season- 
ally adjusted daily water height Danube 
river data. 



ola h P{h) against (^J. 

As it may seen in figure [7] the reverse Log-BHP falls on top of the Log- 
histogram of the deseasonalised Danube data. 

The analogue histograms for the Douro river show a higher concentration of 
values below the origin (mean for the original data). This is not surprising 
because for the Douro basin there is no melting ice to feed the river in 
the spring and in the long dry summer where the natural flow is much 
less than that of winter. The mean daily flow of Douro is close to the 
percentile 70. For the Danube data the average height is about the same as 
the median. The predominance of small runoffs for the Douro may be seen 
in the deseasonalised histogram and Log histogram, figures l8|9l In this last 
graphic the differences for the Danube are made very clear because of the 
spike at the origin. 

Finally, in figure [10] is shown the Log of the BHP and the Gaussian fitd on 
top of the Log histogram of the deseasonalised daily runoff data for the full 
year and only for the winter regime. The BHP PDF comes from a system 



3 only a subset of data was used for the fit. 
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Figure 7: Logarithm of the reversed BHP PDF on top of the Log histogram of the 
Danube seasonally adjusted data. Near the origin there is less deviation from the BHP. 
The left tail of the Log histogram seems to follow a different asymptote than that of the 
BHP. The right tail the of the BHP asymptote seems to be a good fit and the higher 
dispersion around the BHP curve may be due to measurement error because the values 
used in that section are the height maximums. 




Figure 8: Histogram of the seasonally 
adjusted daily water height Douro river 
data. 



Figure 9: Log histogram of the season- 
ally adjusted daily water runoff of the 
Douro river data. The spike at the origin 
(mean of the original data) is explained 
by the high frequency of small runoffs. 
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Figure 10: Log Histogram of the sea- 
sonally adjusted Douro river daily wa- 
ter level, for the full year with the BHP 
and the Gaussian fit on top. The devi- 
ations from the BHP at the origin are 
well fitted with the Gaussian PDF. This 
fact may be due to regulation measures 
at the dams. For smaller and specially 
for larger runoffs the seasonally adjusted 
flow tend to the BHP PDF. 




-2-1 1 2 3 
dh 



Figure 11: Log Histogram of the sea- 
sonally adjusted Douro river daily water 
level, autumn and winter with the BHP 
and the Gaussian fit on top. In these 
seasons the smaller runoffs can be bet- 
ter regulated towards the Gaussian PDF 
and larger runoffs tend to occur more fre- 
quently which explains the larger num- 
ber of points around the BHP curve. 
The resulting PDF is a mixture of a 
gaussian up to a certain control thresh- 
old and the BHP for larger runoffs. 



with no outside tuning so it should model data that is very close to the 
natural flow of the river. The river flow is more close to its natural form 
when regulation is not applicable. This is the case when the river carries 
so much water that the dams have to be set open, see figures [10] and [TTJ 
The differences in the left branch of figures [10] and [11] appears to be caused 
by a regulation for small Winter runoffs which seems not to be possible or 
desirable for the summer resulting in a bad fit for the full year. The PDF of 
the deseasonalised daily Douro runoff data in Winter results is a mixture of 
a Gaussian PDF for small and medium runoffs and BHP for runoffs above 
a certain threshold. Future work passes by searching for stronger evidence 
in favour or against the allegations concerning the differences between the 
so called universal rivers and those which seem not to be universal and also 
the differences between the winter regimes. 

The BHP PDF has been found to be an explanatory model for fluctuations of 
several phenomenon such as, width power in steady state systems, variations 
in river heights and flow, numerous self-organised critical systems among 
others, (see references in |11|). Still it is not known yet why BHP captures 
the essential behaviour of systems of different universality classes. Finally 
there could be a relation between the Fisher-Tippett-Gumbel distribution of 
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extremal statistics and the BHP but recent works [5J 0] , showed asymptotic 
differences between the two so in this direction there is also work to be done. 
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